Numerical method for simulating subsonic flows based on euler equations in lagrangian formulation

ABSTRACT

Numerical methods for simulating subsonic flows and solving inverse problems based on the new two-dimensional Euler equations in Lagrangian formulation. A transformation to derive the Euler equations in Lagrangian plane to simplify the computational grid and to minimize the numerical diffusion caused by the convective flux. Constructed is a numerical scheme, the dimensional-splitting with hybrid upwind flux operators, where the two-dimensional Riemann problem was solved in the Lagrangian plane. The method, solving the inverse geometry shape design problem in Lagrangian plane, gives the flow filed results and the geometry shape concurrently.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a continuation of International Patent ApplicationNo. PCT/CN2010/076768, with an international filing date of Sep. 9,2010, designating the United States, now pending. The contents of all ofthe aforementioned applications, including any intervening amendmentsthereto, are incorporated herein by reference.

BACKGROUND OF THE INVENTION

1. Field of the Invention

This invention relates to the numerical methods for simulating invicidsubsonic flows and solving a class of inverse problems, which is theshape design of aerodynamic body in invicid subsonic flows.

2. Description of the Related Art

Most of the existing CFD (Computational Fluid Dynamics) works utilizethe Eulerian formulation to solve the fluid flow governing equations,which means in the Eulerian plane constituted by the Cartesiancoordinates the computational grid is generated in advance and mainlybased on the geometry constraints. The grid forms the computationalcells. Across the cell interfaces, the convective flux exists since thefluid particles pass over. It is this convective flux that causes severenumerical diffusion in the numerical solutions because the numericaldiffusion is directly associated with the error resulting fromnumerically approximating the convective term. Since the last century,the primary CFD efforts of algorithm researchers had been concentratedon developing more robust, accurate, and efficient ways to reduce thediffusion. In particular, the upwind methods had a great deal of successin solving fluid flows, because the upwind methods reasonably representthe characteristics of the convective flux. Typically, the Godunovmethod^([1]), solving the Riemann problem on the cell interfaces, givesthe most accurate results. The FVS (Flux Vector Splitting) method^([2]),applying the eigenstructure of the equations to treat the flux term onthe cell interfaces, is more efficient. However, the convectivediffusion, as a result of the fundamental nature of the Eulerianformulation, still exists in those characteristic-based numericalschemes.

On the other hand, Lagrangian description to fluid flow states themotions and properties of the given fluid particles as they travel todifferent locations. For the governing equations of the fluid flow, suchas the Euler equations in Lagrangian formulation, there must be onecoordinate representing the streamlines, which can be the streamfunction numerically. Another coordinate could be the distance of theparticles travel. This coordinate system constitutes the Lagrangianplane, where the computational grid points are literally fluid particlesand the computational grid is always simply rectangular. In particular,since the particle paths coincide with the streamlines in steady flowand no fluid particles will cross the streamlines, there is noconvective flux across the cell interfaces and the numerical diffusioncan be thus minimized when the equations are solved in the Lagrangianplane.

One advantage from applying the Lagrangian formulation is when solving aclass of inverse problems. The typical inverse problem of aerodynamicsis posed by specifying the pressure distribution on the solid-wall ofaerodynamic body and determining the geometry of this solid-wall thatrealizes this pressure distribution. If we solve the problems in theEulerian plane, using such as the adjoint method^([3]), we have tofirstly estimate the unspecified body shape, then to generate a gridaround this shape, and to apply the numerical methods to the flows andfind the pressure distribution on the body. Next, a very important andtime-consuming step is to solve the adjoint equation to modify thegeometry shape. As this process has to be repeated for severalsuccessive steps for shape modification until the final target geometryis reached, the computation time of this process is usually long. ALagrangian formulation, as mentioned before, using the stream functionas one coordinate, is more suitable to work on the unspecified geometryproblem, since the unspecified shape (solid-wall boundaries) are alwaysrepresented by streamlines, there no need to modify the computationalgrid like that in the adjoint method. In the Lagrangian plane, thecomputational grid will keep on the same in anytime, no matter thegeometry of the body shape how to change, because any solid-wall surfacein Lagrangian plane is a straight line with the constant values ofstream function. The method for solving the geometry shape designproblem in Lagrangian plane comes to the optimal process (the mostefficient process).

Although the Lagrangian formulation presents such excellent properties,it was limited in supersonic flows simulation before^([4]). When theEuler equations are numerically solved in Lagrangian plane, theyspatially evolve downstream, there is no any upstream information needto consider, which therefore perfectly follows the physicalcharacteristics in supersonic flows. Recently the Lagrangian formulationhad been applied on the aerodynamic shape design problems intwo-dimensional supersonic flow field^([5]). Until recently, for thesubsonic or low speed flow domain, using the advantage of the streamfunction as one coordinate working on the inverse shape design problemsis only for potential flows (invicid, irrotational, incompressible) andlinearized compressible flows^([6]).

There is an apparent need to apply the advantages of the Lagrangianformulation, such as the minimum numerical diffusion in simulation andthe optimal process in shape design. Some industrial applications, suchas two-dimensional invicid subsonic flow simulations and design casesfor nozzle and airfoil, are common, useful and necessary in thepreliminary phase of product design. Numerically solving subsonic flowsby using the strict Lagrangian concept becomes an excessive obstacle,physically because there exist the upstream-propagating waves insubsonic flows. Thus, the existence of a body located downstream istransmitted to the oncoming fluid particles via the waves so that theparticles, sensing the influences of upstream-propagating, can changemotion accordingly. A key to the success of a numerical Lagrangianmethod for subsonic flows lies in how to properly and instantaneouslyfeed these upstream-propagating waves to the particles. Several objectsof this invention are: (a) to provide a Lagrangian formulation of thetwo-dimensional Euler equations with the zeroed-out convective fluxalong one coordinate to minimize numerical diffusion; (b) to provide thenumerical methods that solve the Euler equations in Lagrangianformulation to find more accurate solutions; and (c) to provide thenumerical method that is optimal to solve the inverse shape designproblems in subsonic flows.

SUMMARY OF THE INVENTION

This invention supplies a new technology to simulate the invicidsubsonic flows. It comprises a transformation to transfer thetwo-dimensional Euler equations into the ones in Lagrangian formulationin Lagrangian plane, which will not only accept a simple rectangulargrid, but also produce the minimum convective diffusion when thecharacteristic-based numerical methods are use to solve the equations. Atwo dimensional Riemann solver is created to solve the Riemann problemacross streamlines in Lagrangian plane. An advantage of using this newtechnology lies on solving a class of inverse problem of aerodynamicshape design. Solving the inverse Riemann problem on the solid-wallboundary greatly reduces the complicity of this kind problem, and thecomputing time reduces as much as the order of one. It gives the flowfield solution and the geometry shape concurrently. So solving thisclass of inverse problem comes to the most efficient.

A method for simulating the two-dimensional invicid subsonic flows,comprising the following steps: (i) transforming the two-dimensionalunsteady Euler equations in Euler plane to the Euler equations inLagrangian formulation in Lagrangian plane; (ii) reading object geometryfor providing points on a surface of an object; (iii) establishing acomputational mesh around the object; and (iv) solving the Eulerequations in Lagrangian formulation in Lagrangian plane using theStrang-dimensional-splitting scheme with hybrid upwind flux operators.

Euler equations in Lagrangian formulation include two set of equations,the Lagrangian physical conservation equations and the geometricconservation equations.

The Lagrangian plane has one time coordinate direction, which is statedby the Lagrangian time τ, and two space coordinates directions, which isstated by the stream function ξ and the Lagrangian-distance λ.

The transforming the two-dimensional unsteady Euler equations comprisesa transformation matrix, which is be written as

$J = \begin{bmatrix}1 & 0 & 0 \\{h\; u} & {\cos\;\theta} & U \\{h\; v} & {\sin\;\theta} & V\end{bmatrix}$where h is called Lagrangian behavior parameter, θ is the flow directionangle, u, v are the Cartesian components of flow velocity V, and U, Vare the Lagrangian geometry state variables.

The Lagrangian physical conservation equations read as

${{\frac{\partial f_{L}}{\partial\tau} + \frac{\partial F_{L}}{\partial\lambda} + \frac{\partial G_{L}}{\partial\xi}} = 0},$where f_(L) is the conservation variables vector and F_(L), G_(L) arethe convective fluxes respectively in the λ, and ξ direction inLagrangian plane,

${f_{L} = \begin{bmatrix}{\rho\; J} \\{\rho\; J\; u} \\{\rho\; J\; v} \\{\rho\; J\; E}\end{bmatrix}},\mspace{14mu}{F_{L} = \begin{bmatrix}{\left( {1 - h} \right)\rho\; K} \\{{\left( {1 - h} \right)\rho\; K\; u} + {V\; p}} \\{{\left( {1 - h} \right)\rho\; K\; v} - {U\; p}} \\{\left( {1 - h} \right)\rho\; K\; H}\end{bmatrix}},{G_{L} = \begin{bmatrix}0 \\{{- p}\;\sin\;\theta} \\{p\;\cos\;\theta} \\0\end{bmatrix}},\mspace{14mu}{\theta = {t\;{g^{- 1}\left( \frac{v}{u} \right)}}},\mspace{14mu}{K = {{u\; V} - {v\; U}}},\mspace{14mu}{0 < h \leq 1.}$where the variables t, ρ, p, E and H are respectively the time, fluiddensity, pressure, total energy and enthalpy.

The geometric conservation equations read as

${{\frac{\partial f_{g}}{\partial\tau} + \frac{\partial F_{g}}{\partial\lambda} + \frac{\partial G_{g}}{\partial\xi}} = 0},$where

${f_{g} = \begin{bmatrix}U \\V\end{bmatrix}},\mspace{14mu}{F_{g} = \begin{bmatrix}0 \\0\end{bmatrix}},\mspace{14mu}{G_{g} = {\begin{bmatrix}{{- h}\; u} \\{{- h}\; v}\end{bmatrix}.}}$

The computing mesh is the two-dimensional rectangular mesh with λ and ξas the two coordinate directions in the Lagrangian plane.

Solving the Euler equations in Lagrangian formulation in Lagrangianplane utilizes the time-like iteration along the direction of theLagrangian time τ to find the steady solutions of the equations.

The Strang-dimensional-splitting scheme with hybrid upwind fluxoperators comprises the following steps in one step of updating of theconservative variables when the Lagrangian time τ in Lagrangian planeevolves, (A) calculating the convective fluxes across the cellinterfaces in λ-direction by the FVS method with the time-step Δτ; (B)solving the Riemann problem across streamlines in ξ-direction with thehalf time-step Δτ by the Godunov method with the local Riemann solversequentially; and (C) calculating the convective fluxes across the cellinterfaces in λ-direction by the FVS method with the time-step Δτ again.

The Riemann problem across streamlines in ξ-direction comprises that atthe left or right state (region) of the cell interface, could be shockand expansion waves respectively, between which it is the star state(region), which can be divided into left and right star states (region)again.

Solving the Riemann problem across streamlines comprise the followingsteps: 1) Connecting the left and right states to the star state byintegrating along the characteristic equations of the Lagrangianphysical conservation equations; 2) Recovering the velocity magnitude inthe star state; 3) Solving the combination function to find the flowangle in the star state; and 4) Finding the velocity component in thestar state.

The step of recovering the velocity magnitude in the star region can bedone according to the Rankine-Hugoniot relations across shocks and theEnthalpy constants across expansion waves.

The combination function f(u, v) reads as

${f\left( {u,v} \right)} = {{\frac{1}{2}\frac{v\sqrt{u^{2} + v^{2}}}{u}} + {\frac{1}{2}u\;{{\ln\left( {v + \sqrt{u^{2} + v^{2}}} \right)}.}}}$

A method for solving the inverse shape design problem, comprising thefollowing steps: a) initializing the flow field parameters andsolid-wall angle; b) inputting the specified pressure distribution onthe solid-wall; c) solving the inversed Riemann problem to find themirror state of the solid-wall boundary; d) calculating the solid-wallangle; e) checking the convergence of the solid-wall angle; f) updatingthe flow field parameters; and g) repeating (c).

Solving the inversed Riemann problem comprises solving the combinationfunction based on the specified pressure on the solid-wall, which can bewritten as

${{f\left( {V_{R},\theta_{L}} \right)} = {{f\left( {V_{R},\theta_{R}} \right)} + \frac{2\left( {P_{w} - P_{R}} \right)}{\rho_{R}a_{R}}}},$where P_(R), ρ_(R), a_(R), V_(R), θ_(R), are the known flow fieldparameters in the solid-wall boundary, f is the combination function,p_(w) is the specified pressure, and θ_(L) is the flow angle in ghostcells.

Solving the inversed Riemann problem comprises solving the combinationfunction based on the specified pressure on the solid-wall by theNewton-Raphson iteration method to find the flow angle θ_(L).

BRIEF DESCRIPTION OF THE DRAWINGS

A detailed description of accompanying drawings will be provided below:

FIG. 1 illustrates the Riemann problem across streamline alongξ-direction in Lagrangian plane;

FIG. 2 illustrates flow chart of one time-step for updating the flowfield using the second-order Scheme: dimensional-splitting with hybridflux operators;

FIG. 3 illustrates flow chart of an inverse method for geometry shapedesign;

FIG. 4 illustrates the Lagrangian and Eulerian solutions for parabolicnozzle flow at M_(in)=0.5

-   -   (a) Grid in Lagrangian plane,    -   (b) Grid in Eulerian plane,    -   (c) Streamlines by the Lagrangian solution,    -   (d) Streamlines by the Eulerian solution,    -   (e) Comparison of the pressure distributions by Lagrangian and        exact solution on the bottom solid-wall boundary, and    -   (f) Grid built by the Lagrangian solution; and

FIG. 5 calculated geometry shape of the solid-wall based on theprescribed pressure distribution

-   -   (a) Grid in Eulerian plane,    -   (b) Pressure distribution on the solid wall, and    -   (c) Designed solid-wall.

REFERENCE CHARACTERS

$\begin{matrix}{{1.\;\lbrack Q\rbrack}_{i,j}^{n},} & \; \\{{2.\;\lbrack Q\rbrack}_{{i \pm {1/2}},j}^{n},} & \; \\{{3.\;\lbrack F\rbrack}_{{i \pm {1/2}},j}^{n},} & \; \\{{{4.\mspace{11mu} f_{i,j}^{n + {1/3}}} = {f_{i,j}^{n} - {{\Delta\tau}\left\lbrack \frac{\left( {F_{{i + {1/2}},j}^{n} - F_{{i - {1/2}},j}^{n}} \right)}{\Delta\lambda} \right\rbrack}}},} & \; \\{{5.\mspace{11mu}(Q)_{i,j}^{n + {1/3}}},} & \; \\{{6.\;\lbrack Q\rbrack}_{i,{j \pm {1/2}}}^{n + {1/3}},} & \; \\{{7.\;\lbrack G\rbrack}_{{i \pm {1/2}},j}^{n + {1/3}},} & \; \\{{{8.\mspace{11mu} f_{i,j}^{n + {2/3}}} = {f_{i,j}^{n + {1/3}} - {{\Delta\tau}\left\lbrack \frac{\left( {G_{i,{j + {1/2}}}^{n + {1/3}} - G_{i,{j - {1/2}}}^{n + {1/3}}} \right)}{\Delta\xi} \right\rbrack}}},} & \; \\{{9.\;\lbrack Q\rbrack}_{i,j}^{n + {2/3}},} & \; \\{{10.\;\lbrack Q\rbrack}_{{i \pm {1/2}},j}^{n + {2/3}},} & \; \\{{11.\;\lbrack F\rbrack}_{{i \pm {1/2}},j}^{n},} & \; \\{{{12.\mspace{11mu} f_{i,j}^{n + 1}} = {f_{i,j}^{n + {2/3}} - {{\Delta\tau}\left\lbrack \frac{\left( {F_{{i + {1/2}},j}^{n + {2/3}} - F_{{i - {1/2}},j}^{n + {2/3}}} \right)}{\Delta\lambda} \right\rbrack}}},} & \; \\{{13.\;\lbrack Q\rbrack}_{i,j}^{n + 1}.} & \;\end{matrix}$

DETAILED DESCRIPTION OF THE EMBODIMENTS

The Two-Dimensional Euler Equations in Lagrangian Formulation

This invention firstly provides a coordinate transformation from theEulerian variables (t, x, y) in the Eulerian plane to the Lagrangianvariables (r, λ, ξ) in the Lagrangian plane, where the variable τ withthe same dimension as time term t in the Euler plane, is introduced tofunction as the ‘time-like’ marching term for subsonic flows; the othertwo independent variables are the stream function ξ with dimension[m²·s⁻¹] and the Lagrangian-distance λ with dimension [m]. Thecoordinate ξ shall coincide with the streamlines of the fluid particlesand λ is specified as the location of different particles alongstreamlines. This invention stars from the Euler equations indifferential conservation form for two-dimensional unsteady inviscidcompressible flows in the Eulerian plane

$\begin{matrix}{{{\frac{\partial f}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y}} = 0},} & (1)\end{matrix}$Where f is the conservation variables vector and F, G are the convectivefluxes in the two Cartesian directions. In detail,

$\begin{matrix}{{f = \begin{bmatrix}\rho \\{\rho\; u} \\{\rho\; v} \\{\rho\; E}\end{bmatrix}},\mspace{14mu}{F = \begin{bmatrix}{\rho\; u} \\{{\rho\; u^{2}} + \rho} \\{\rho\; u\; v} \\{\left( {{\rho\; E} + p} \right)u}\end{bmatrix}},\mspace{14mu}{G = {\begin{bmatrix}{\rho\; v} \\{\rho\; u\; v} \\{{\rho\; v^{2}} + p} \\{\left( {{\rho\; E} + p} \right)v}\end{bmatrix}.}}} & (2)\end{matrix}$

The variables t, u, v, ρ and p are respectively the time, Cartesiancomponents of flow velocity V, fluid density and pressure. The totalenergy E and enthalpy H are

$\begin{matrix}{{E = {{\frac{1}{2}\left( {u^{2} + v^{2}} \right)} + {\frac{1}{r - 1}\frac{p}{\rho}\mspace{14mu}{and}}}}\mspace{14mu}{H = {{E + \frac{p}{\rho}} = {\frac{u^{2} + v^{2}}{2} + {\frac{r}{r - 1}{\frac{p}{\rho}.}}}}}} & (3)\end{matrix}$

In above equation, γ is the specific heat ratio. The Jacobian matrix ofthe transformation from coordinates (t, x, y) to (τ, λ, ξ) can bewritten as

$\begin{matrix}{{J = {\frac{\partial\left( {t,x,y} \right)}{\partial\left( {\tau,\lambda,\xi} \right)} = \begin{bmatrix}1 & 0 & 0 \\{h\; u} & {\cos\;\theta} & U \\{h\; v} & {\sin\;\theta} & V\end{bmatrix}}},} & (4)\end{matrix}$and its determinant J becomes

$\begin{matrix}{{J = {{\begin{matrix}{\cos\;\theta} & U \\{\sin\;\theta} & V\end{matrix}} = {{V\;\cos\;\theta} - {U\;\sin\;\theta}}}},} & (5)\end{matrix}$where h is called Lagrangian behavior parameter, θ is the flow directionangle and the two quantities called the Lagrangian geometry statevariables U, V with the dimension [m²s·kg⁻¹] are defined as

$\begin{matrix}{{U = \frac{\partial x}{\partial\xi}},\mspace{14mu}{V = {\frac{\partial y}{\partial\xi}.}}} & (6)\end{matrix}$

Finally, the two-dimensional time-dependent Euler equations (1) havebeen transformed into the ones in Lagrangian formulation in theLagrangian plane, which include two set of partial differentialequations in the conservation form,

$\begin{matrix}{{{\frac{\partial f_{L}}{\partial\tau} + \frac{\partial F_{L}}{\partial\lambda} + \frac{\partial G_{L}}{\partial\xi}} = 0},} & (7) \\{{{\frac{\partial f_{g}}{\partial\tau} + \frac{\partial F_{g}}{\partial\lambda} + \frac{\partial G_{g}}{\partial\xi}} = 0},} & (8)\end{matrix}$where f_(L), F_(L) and G_(L) have the same definition as f, F, G in (2)but in Lagrangian plane,

$\begin{matrix}{{f_{L} = \begin{bmatrix}{\rho\; J} \\{\rho\; J\; u} \\{\rho\; J\; v} \\{\rho\; J\; E}\end{bmatrix}},\mspace{14mu}{F_{L} = \begin{bmatrix}{\left( {1 - h} \right)\rho\; K} \\{{\left( {1 - h} \right)\rho\; K\; u} + {V\; p}} \\{{\left( {1 - h} \right)\rho\; K\; v} - {U\; p}} \\{\left( {1 - h} \right)\rho\; K\; H}\end{bmatrix}},\mspace{14mu}{G_{L} = \begin{bmatrix}0 \\{{- p}\;\sin\;\theta} \\{p\;\cos\;\theta} \\0\end{bmatrix}},} & (9) \\{{f_{g} = \begin{bmatrix}U \\V\end{bmatrix}},\mspace{14mu}{F_{g} = \begin{bmatrix}0 \\0\end{bmatrix}},\mspace{14mu}{G_{g} = \begin{bmatrix}{{- h}\; u} \\{{- h}\; v}\end{bmatrix}},} & (10)\end{matrix}$with the notations of J in (5), and

$\begin{matrix}{{\theta = {{tg}^{- 1}\left( \frac{v}{u} \right)}},} & (11) \\{{K = {{uV} - {vU}}},} & (12) \\{0 < h \leq 1.} & (13)\end{matrix}$

The equations (7) are called the Lagrangian physical conservationequations and (8) are the geometric conservation equations. In (9), wefind that the first and last element in vector G_(L) are zero, whichmeans that for the computing grid in Lagrangian plane there apparentlyare no flux across the cell interfaces along ξ-direction for thecontinuity and energy equations, so that the numerical diffusion isminimized. There are some other properties of the Euler equations inLagrangian formulation.

-   -   1. The equations (7) satisfy the homogeneity property, which        implies that the FVS method can be used to resolve the flux        term.    -   2. The equations (7) are not with the rotational invariance        property, which implies that it is better to use the different        numerical schemes for the convective flux along the two        coordinate directions to keep the characteristic non-smeared.    -   3. The partial differential equation systems (7) and (8) are        both strongly-hyperbolic, which implies that it is the        well-posed initial problem and they can be solved using the        time-like marching scheme.

Eigenstructure and Characteristics of the Euler Equations in LagrangianFormulation

For the physical conservation equations (7), the eigenstructure (eigenvalues and eigen vectors) can be found, according to the definition ofthe Jacobian matrix,

$\begin{matrix}{{B_{L} = \frac{\partial F_{L}}{\partial Q}},\mspace{14mu}{C_{L} = \frac{\partial G_{L}}{\partial Q}},} & (14)\end{matrix}$where Q=[ρ, u, v, p]^(T). Ones can find the four eigenvalues of theJacobian matrix B_(L),μ_(L) _(1,2) =(1−h)√{square root over (u ² +v ²)}μ_(L) _(3,4)=(1−h)√{square root over (u ² +v ²)}±(a/j)√{square root over (U ² +V²)}  (15)and the four corresponding right eigenvectors

$\begin{matrix}{{K_{L}^{(1)} = \begin{bmatrix}1 \\0 \\0 \\0\end{bmatrix}},{K_{L}^{(2)} = \begin{bmatrix}0 \\{U/\rho} \\{V/\rho} \\0\end{bmatrix}},{K_{L}^{(3)} = \begin{bmatrix}1 \\{\frac{V}{\sqrt{U^{2} + V^{2}}}\left( \frac{a}{\rho} \right)} \\{\frac{- U}{\sqrt{U^{2} + V^{2}}}\left( \frac{a}{\rho} \right)} \\a^{2}\end{bmatrix}},{K_{L}^{(4)} = \begin{bmatrix}1 \\{\frac{- V}{\sqrt{U^{2} + V^{2}}}\left( \frac{a}{\rho} \right)} \\{\frac{U}{\sqrt{U^{2} + V^{2}}}\left( \frac{a}{\rho} \right)} \\a^{2}\end{bmatrix}}} & (16)\end{matrix}$along λ-direction; and the four eigenvalues of the Jacobian matrix C_(L)μ_(L) _(1,2) =0, μ_(L) _(3,4) =±(a/j),  (17)and the four corresponding right eigenvectors

$\begin{matrix}{{K_{L}^{(1)} = \begin{bmatrix}1 \\0 \\0 \\0\end{bmatrix}},{K_{L}^{(2)} = \begin{bmatrix}0 \\{\cos\;{\theta/\rho}} \\{\sin\;{\theta/\rho}} \\0\end{bmatrix}},{K_{L}^{(3)} = \begin{bmatrix}1 \\{{- \sin}\;{\theta\left( \frac{a}{\rho} \right)}} \\{\cos\;{\theta\left( \frac{a}{\rho} \right)}} \\a^{2}\end{bmatrix}},{K_{L}^{(4)} = \begin{bmatrix}1 \\{\sin\;{\theta\left( \frac{a}{\rho} \right)}} \\{{- \cos}\;{\theta\left( \frac{a}{\rho} \right)}} \\a^{2}\end{bmatrix}},} & (18)\end{matrix}$in ξ-direction, where a is the speed of sound. The characteristicequations of the Lagrangian physical conservation equations (7) can bewritten as

$\begin{matrix}{{{{dp} \pm {\rho\; a\frac{\sqrt{U^{2} + V^{2}}}{V}{du}}} = 0}{along}\;{{\frac{\mathbb{d}\lambda}{\mathbb{d}\tau} = {\mu_{L_{3,4}} = {{\left( {1 - h} \right)\sqrt{u^{2} + v^{2}}} \pm {\left( \frac{a}{J} \right)\sqrt{{U^{2} + V^{2}}\;}}}}};}} & (19) \\{{{dp} \pm {\rho\; a\frac{\sqrt{u^{2} + v^{2}}}{u}{dv}}} = 0} & (20)\end{matrix}$along dξ/dτ=μ_(L) _(3,4) =±(a/j). Based on the Jacobian matrix of thetransformation (4), the Riemann invariants in Lagrangian plane inλ-direction become

u - V U 2 + V 2 ⁢ ⁢ 2 ⁢ ⁢ a γ - 1 = + , ⁢ u + V U 2 + V 2 ⁢ ⁢ 2 ⁢ ⁢ a γ - 1 = - ;( 21 )the Riemann invariants in ξ-direction are

v - cos ⁢ ⁢ θ ⁢ 2 ⁢ ⁢ a γ - 1 = + , ⁢ v + cos ⁢ ⁢ θ ⁢ 2 ⁢ ⁢ a γ - 1 = - . ( 22 )

The equations (15)-(22) will be used to construct the numerical schemefor solving the governing equations (7) and (8).

Construction of Numerical Scheme

To solve the two-dimensional Euler equations in Lagrangian formulation(7) and (8) in Lagrangian plane, the FVM (Finite Volume Method) ischosen as the spatial discretization. As stated before, thecomputational grid in Lagrangian plane is simple rectangular. Acell-center scheme will be applied, where the flow quantities are storedat the centers of the grid cells. Based on the properties of thehomogeneity, rotation invariance and hyperbolicity of the equations (7)and (8), to keep the efficiency and accuracy of the calculations, thisinvention, combining the advantages of both the Godunov method and FVSmethod, supplies a new method of solution, the dimensional splittingscheme with hybrid upwind flux operators, which means that theconvective fluxes across the cell interfaces will be calculated by theFVS method in λ-direction and the Godunov method with local Riemannsolver in ξ-direction sequentially. The Strang dimensional-splitting,which is the second-order accuracy in time, is applied for the temporaldiscretization. Now, let's drop the subscript L with the Lagrangiansense for the convenience reason.

FVS Method in λ-Direction

Firstly rewrite equations (7) only along λ-direction,

$\begin{matrix}{{\frac{\partial f}{\partial\tau} + \frac{\partial F}{\partial\lambda}} = 0.} & (23)\end{matrix}$

By applying the FVS (flux vector splitting) method, one obtains thesolutions at all the cell (i,j) at the new time-step

$\begin{matrix}{{f_{i,j}^{n + 1} = {f_{i,j}^{n} - {\frac{\Delta\tau}{\Delta\lambda}\left\lbrack {{L_{i}^{+}f_{i,j}^{n}} - {L_{i - 1}^{+}f_{{i - 1},j}^{n}} + {L_{i + 1}^{-}f_{{i + 1},j}^{n}} - {L_{i}^{-}f_{i,j}^{n}}} \right\rbrack}}},} & (24)\end{matrix}$where n is the number of time-step, i,j are the cell indices, Δτ thetime-step, Δλ is the spatial step. The split matrixL ⁺ =K _(Lcv)Λ⁺ K _(Lcv) ⁻¹ and L ⁻ =K _(Lcv)Λ⁻ K _(Lcv) ⁻¹,  (25)with the split matrix of eigenvalues

$\begin{matrix}{{\Lambda^{+} = {{\begin{bmatrix}\mu_{L_{1}} & \; & \; & \; \\\; & \mu_{L_{2}} & \; & \; \\\; & \; & \mu_{L_{3}} & \; \\\; & \; & \; & 0\end{bmatrix}\mspace{14mu}{and}\mspace{14mu}\Lambda^{-}} = \begin{bmatrix}0 & \; & \; & \; \\\; & 0 & \; & \; \\\; & \; & 0 & \; \\\; & \; & \; & \mu_{L_{4}}\end{bmatrix}}};} & (26)\end{matrix}$and the matrix

$\begin{matrix}{K_{Lcv} = \begin{bmatrix}1 & 0 & 1 & 1 \\u & {U\frac{\rho}{J}} & {u + \frac{Va}{U^{2} + V^{2}}} & {u - \frac{Va}{U^{2} + V^{2}}} \\v & {V\frac{\rho}{J}} & {v - \frac{Ua}{U^{2} + V^{2}}} & {v + \frac{Ua}{U^{2} + V^{2}}} \\\frac{u^{2} + v^{2}}{2} & {\left( {{Uu} + {Vv}} \right)\frac{\rho}{J}} & {\frac{u^{2} + v^{2}}{2} + \frac{Ka}{\sqrt{U^{2} + V^{2}}} + \frac{a^{2}}{1 - \gamma}} & {\frac{u^{2} + v^{2}}{2} - \frac{Ka}{\sqrt{U^{2} + V^{2}}} + \frac{a^{2}}{1 - \gamma}}\end{bmatrix}} & (27)\end{matrix}$and the inverse of it reads

$\begin{matrix}{K_{Lcv}^{- 1} = {{\frac{\left( {\gamma - 1} \right)}{2\; a^{2}}\begin{bmatrix}{u^{2} + v^{2} + \frac{2\; a^{a}}{\gamma - 1}} & {{- 2}\; u} & {{- 2}\; v} & 2 \\{{- \left( \frac{{Uu} + {Vv}}{U^{2} + V^{2}} \right)}\frac{2\; a^{2}J}{\left( {\gamma - 1} \right)\rho}} & {\left( \frac{U}{U^{2} + V^{2}} \right)\frac{2\; a^{2}J}{\left( {\gamma - 1} \right)\rho}} & {\left( \frac{V}{U^{2} + V^{2}} \right)\frac{2\; a^{2}J}{\left( {\gamma - 1} \right)\rho}} & 0 \\{{- \frac{u^{2} + v^{2}}{2}} - {\frac{K}{\sqrt{U^{2} + V^{2\;}}}\left( \frac{a}{\gamma - 1} \right)}} & {u + {\frac{V}{\sqrt{U^{2} + V^{2\;}}}\left( \frac{a}{\gamma - 1} \right)}} & {v - {\frac{U}{\sqrt{U^{2} + V^{2\;}}}\left( \frac{a}{\gamma - 1} \right)}} & {- 1} \\{{- \frac{u^{2} + v^{2}}{2}} + {\frac{K}{\sqrt{U^{2} + V^{2\;}}}\left( \frac{a}{\gamma - 1} \right)}} & {u - {\frac{V}{\sqrt{U^{2} + V^{2\;}}}\left( \frac{a}{\gamma - 1} \right)}} & {v + {\frac{U}{\sqrt{U^{2} + V^{2\;}}}\left( \frac{a}{\gamma - 1} \right)}} & {- 1}\end{bmatrix}}.}} & (28)\end{matrix}$

Godunov Method with Riemann Solver in ξ-Direction

For the Godunov method, the key ingredient is the solution of theRiemann problem.

This invention solves the Riemann problem along ξ-direction, which iscalled the Riemann problem across streamlines. FIG. 1 shows thedefinition of the Riemann problem across streamlines. Rewrite (7) onlyin ξ-direction,

$\begin{matrix}{{{\frac{\partial f}{\partial\tau} + \frac{\partial{G(f)}}{\partial\xi}} = 0},\mspace{11mu}{f = \left\{ \begin{matrix}f_{L} & {\xi < 0} \\f_{R} & {\xi > 0.}\end{matrix} \right.}} & (29)\end{matrix}$

Where f_(L) and f_(R) are the conservation variables at the left andright state of the cell interface respectively. From now, the subscriptsL are for the left state, R for the right, and * for the star for theRiemann problem. At the left or right state, shock and expansion wavescould be. Between them, it is the star state (region), which can bedivided into left and right star states again like shown in FIG. 1. Thesolutions of the new time-step can be obtained from the discretizedequations of (29),

$\begin{matrix}{f_{i,j}^{n + \frac{1}{2}} = {f_{i,j}^{n} - {{\frac{\Delta\tau}{\Delta\xi}\left\lbrack {{G\left( f_{i,{j + \frac{1}{2}}}^{n} \right)} - {G\left( f_{i,{j + \frac{1}{2}}}^{n} \right)}} \right\rbrack}.}}} & (30)\end{matrix}$

The task is to find the values of the primitive variables p, ρ, u, vfrom f to get G at cell interfaces j±½, where the primitive variablesare the variables in the star region in the Riemann problem. A Riemannsolver was built as the following procedure. Connecting the left andright states to the star state by integrating along the characteristicsequations (19) and (20), then ones have

$\begin{matrix}\left\{ {{{\begin{matrix}{{{p_{*} + {C_{L} \cdot {f\left( {u_{*},v_{*}} \right)}}} = {p_{L} + {C_{L} \cdot {f\left( {u_{L},v_{L}} \right)}}}},} \\{{{p_{*} - {C_{R} \cdot {f\left( {u_{*},v_{*}} \right)}}} = {p_{R} - {C_{R} \cdot {f\left( {u_{R},v_{R}} \right)}}}},(32)} \\{{{\rho_{*L} + {{f\left( {u_{*},v_{*}} \right)}\left( \frac{\rho_{L}}{a_{L}} \right)}} = {\rho_{L} + {{f\left( {u_{L},v_{L}} \right)}\left( \frac{\rho_{L}}{a_{L}} \right)}}},(33)} \\{{{\rho_{*R} - {{f\left( {u_{*},v_{*}} \right)}\left( \frac{\rho_{R}}{a_{R}} \right)}} = {\rho_{R} - {{f\left( {u_{R},v_{R}} \right)}\left( \frac{\rho_{R}}{a_{R}} \right)}}},(34)}\end{matrix}{where}C_{L}} = {\rho_{L}a_{L}}},{C_{R} = {\rho_{R}a_{R}}},{and}} \right. & (31) \\{{f\left( {u,v} \right)} = {{\frac{1}{2}\frac{v\sqrt{u^{2} + v^{2}}}{u}} + {\frac{1}{2}u\;{\ln\left( {v + \sqrt{u^{2} + v^{2}}} \right)}}}} & (35)\end{matrix}$f(u, v) is called the combination function, which can be considered asthe combination of the velocity component (u, v) in Lagrangian plane andhave the same role as the velocity term with dimension [m/s]. The (35)has its polar form as

$\begin{matrix}{{f_{\Delta}\left( {V,\theta} \right)} = {{\frac{V}{2}\left( {{{tg}\;\theta} + {\cos\;{{\theta ln}\left( {1 + {\sin\;\theta}} \right)}}} \right)} + {\frac{\cos\;\theta}{2}V\;\ln\;{V.}}}} & (36)\end{matrix}$

With V=√{square root over (u²+v²)}, u=V cos θ and v=V sin θ. Thesolutions of (31)-(34) present the primitive variables in the star state

$\begin{matrix}\left\{ \begin{matrix}{{p_{*} = \frac{{C_{R}p_{L}} + {C_{L}p_{R}} + {C_{L}{C_{R}\left\lbrack {{f\left( {u_{L},v_{L}} \right)} - {f\left( {u_{R},v_{R}} \right)}} \right\rbrack}}}{C_{L} + C_{R}}},} \\{{{f\left( {u_{*},v_{*}} \right)} = \frac{{C_{L} \cdot {f\left( {u_{L},v_{L}} \right)}} + {C_{R} \cdot {f\left( {u_{R},v_{R}} \right)}} + \left\lbrack {p_{L} - p_{R}} \right\rbrack}{C_{L} + C_{R}}},} \\{{\rho_{*L} = {\rho_{L} + {\left( {{f\left( {u_{L},v_{L}} \right)} - {f\left( {u_{*},v_{*}} \right)}} \right)\left( {\rho_{L}/a_{L}} \right)}}},} \\{\rho_{*R} = {\rho_{R} + {\left( {{f\left( {u_{*},v_{*}} \right)} - {f\left( {u_{R},v_{R}} \right)}} \right){\left( {\rho_{R}/a_{R}} \right).}}}}\end{matrix} \right. & \begin{matrix}(37) \\\; \\(38) \\\; \\(39) \\(40)\end{matrix}\end{matrix}$

To find the velocity components (u_(*), v_(*)) either in left or rightstar state, the equation (38) need to be solved. Firstly, one canrewrite (38) asF(u _(*) ,v _(*))=f(u _(*) ,v _(*))−B=0,  (41)with the known conditions in left and right states,

$\begin{matrix}{B = {\frac{{C_{L} \cdot {f\left( {u_{L},v_{L}} \right)}} + {C_{R} \cdot {f\left( {u_{R},v_{R}} \right)}} + \left\lbrack {p_{L} - p_{R}} \right\rbrack}{C_{L} + C_{R}}.}} & (42)\end{matrix}$

The combination function (35) can be further modified. One can definetgθ _(*)=ε,  (43)Where θ_(*) is the flow angle either in left or right star state, then

$\begin{matrix}{{u_{*} = {{V_{*}\cos\;\theta_{*}} = \frac{V_{*}}{\sqrt{1 + ɛ^{2}}}}},\mspace{14mu}{v_{*} = {{V_{*}\sin\;\theta_{*}} = \frac{ɛ\; V_{*}}{\sqrt{1 + ɛ^{2}}}}},} & (44)\end{matrix}$where V_(*) representing the velocity magnitude either in the left orright star state. Further, the equation (41) can be finalized as afunction of ε

$\begin{matrix}{{F(ɛ)} = {{{\frac{V_{*}}{2}\left( {ɛ + \frac{{\ln\; V_{*}} + {\ln\left( {ɛ + \sqrt{1 + ɛ^{2}}} \right)} - {\ln\sqrt{1 + ɛ^{2}}}}{\sqrt{1 + ɛ^{2}}}} \right)} - B} = 0.}} & (45)\end{matrix}$

The work needs to do is numerically to solve the equation F(ε)=0 with agiven velocity V_(*) to find ε for θ_(*). The equation (45) does havethe differentiability and its first-order derivative of F is

$\begin{matrix}{{F^{\prime}(ɛ)} = {\frac{V_{*}}{2}{\left( \frac{\begin{matrix}{{\left( {2 + ɛ^{2}} \right)\sqrt{1 + ɛ^{2}}} - ɛ -} \\{ɛ\left( {{\ln\; V_{*}} + {\ln\left( {ɛ + \sqrt{1 + ɛ^{2}}} \right)} - {\ln\sqrt{1 + ɛ^{2}}}} \right)}\end{matrix}}{\sqrt{\left( {1 + ɛ^{2}} \right)^{3}}} \right).}}} & (46)\end{matrix}$

Numerical experiments show that for the non-dimensional velocity V_(*)≦5(subsonic flows) and ε∈[−1, 1], F′(ε)>0, which means that the functionF(ε) is monotone within this domain; in the mean time, F(−ε)·F(ε)<0, sothat the Newton-Raphson iteration as a root-finding algorithm is a goodway to find the root ε for the equation (45) with a given V_(*). To findthe value of V_(*), the non-linear waves (shear, shock and expansionwaves) in Riemann problem have to be recovered with the approximatevalues of p*, ρ_(*L), ρ_(*R). In this invention, the following relationsare considered.

Rankine-Hugoniot Relations Across Shocks

If a shock wave appears between one state (say left) and thecorresponding star state, across this shock, the Rankine-Hugoniotrelations can be used to the left state and star state. They are,

$\begin{matrix}\left\{ \begin{matrix}{{{\mu_{s}\left( {{\rho_{L}J_{L}} - {\rho_{*L}J_{*L}}} \right)} = 0},} \\{{{\mu_{s}\left( {{\rho_{L}J_{L}u_{L}} - {\rho_{*L}J_{*L}u_{*L}}} \right)} = {{p_{*}\sin\;\theta_{*}} - {p_{L}\sin\;\theta_{L}}}},} \\{{{\mu_{s}\left( {{\rho_{L}J_{L}v_{L}} - {\rho_{*L}J_{*L}v_{*L}}} \right)} = {{p_{L}\cos\;\theta_{L}} - {p_{*}\cos\;\theta_{*}}}},} \\{{{\mu_{s}\left( {{\rho_{L}J_{L}E_{L}} - {\rho_{*L}J_{*L}E_{*L}}} \right)} = 0},}\end{matrix} \right. & \begin{matrix}(46) \\(47) \\(48) \\(49)\end{matrix}\end{matrix}$where μ_(s) is the shock wave speed, J_(*L) and E_(*L), are J(transformation Jacobian) and E (total energy) in the star staterespectively. From (46) and (49), one receives

$\begin{matrix}{{V_{*L} = \sqrt{V_{L} + {\frac{1}{\gamma - 1}\left( {\frac{P_{L}}{\rho_{L}} - \frac{P_{*}}{\rho_{*L}}} \right)}}}{and}} & (50) \\{V_{*R} = {\sqrt{V_{R} + {\frac{1}{\gamma - 1}\left( {\frac{P_{R}}{\rho_{R}} - \frac{P_{*}}{\rho_{*R}}} \right)}}.}} & (51)\end{matrix}$

Velocity absolute value V_(*L) or V_(*R) can be substituted into (45) toreceive θ_(*R) or θ_(*R).

Enthalpy Constants Across Expansion Waves

If an expansion wave appears between one state (say left) and thecorresponding star state, the connection between the two states can befound by considering the enthalpy constant across the waves. Ones have

$\begin{matrix}{{H_{L} = {{{\frac{1}{2}V_{L}} + {\frac{\gamma}{\gamma - 1}\frac{P_{L}}{\rho_{L}}}} = {{{\frac{1}{2}V_{*L}} + {\frac{\gamma}{\gamma - 1}\frac{P_{*}}{\rho_{*L}}}} = H_{*L}}}},} & (52)\end{matrix}$from which the velocity magnitude can be found as

$\begin{matrix}{{V_{*L} = \sqrt{{2H_{L}} - {\frac{2\gamma}{\gamma - 1}\frac{P_{*}}{\rho_{*L}}}}}{and}} & (53) \\{V_{*R} = {\sqrt{{2H_{R}} - {\frac{2\gamma}{\gamma - 1}\frac{P_{*}}{\rho_{*R}}}}.}} & (54)\end{matrix}$

The selection of the non-linear waves is based on the wave-typeconditions, in which the values of pressure in left, right and starstates can be found to judge what non-linear wave could appearing in theRiemann problem. AssumingP _(min)=min(P _(L) ,P _(R));  (55)P _(max)=max(P _(L) ,P _(R)),  (56)which, as well as p_(*), form the following wave-type choosingconditions:

-   -   (1) If P_(min)<P_(*)<P_(max), which means in the Riemann problem        one shock wave in the state with p_(min) and one expansion wave        in the state with max P_(max).    -   (2) If P_(*)≧P_(max), which means there are two shock waves in        Riemann problem.    -   (3) If P_(*)≦P_(min), which means there are two expansion waves        the in Riemann problem.

After θ_(*) is obtained, we can find G_(L) according to (9) and u_(*L),v_(*L) or u_(*R), v_(*R) from (44). Those values, as well as p_(*), canbe used to compute the flux term in Godunov method (30). At the sametime, the geometry conservation equations (8) can be updated

$\begin{matrix}{\begin{bmatrix}U_{i,j}^{n + 1} \\V_{i,j}^{n + 1}\end{bmatrix} = {\begin{bmatrix}U_{i,j}^{n} \\V_{i,j}^{n}\end{bmatrix} + {h\;{{\frac{{\Delta\tau}_{i,j}^{n}}{\Delta\;\xi_{i,j}}\begin{bmatrix}{u_{i,{j + {1/2}}}^{n + {1/2}} - u_{i,{j - {1/2}}}^{n + {1/2}}} \\{v_{i,{j + {1/2}}}^{n + {1/2}} - v_{i,{j - {1/2}}}^{n + {1/2}}}\end{bmatrix}}.}}}} & (57)\end{matrix}$

Although this step is an additional to the method in the Eulerianformulation, it is a trivial step and will no cost more computing time.

Boundary Conditions

Two layers of cells should be put, which are called the ghost cells,surrounding the whole computational domain, to make it possible toextend any spatial discretization scheme beyond the boundaries.

Solid-Wall Boundary Condition

In the selected scheme, the Riemann problem at the solid-wall boundaryin ξ-direction has to be solved, where the left or right state inRiemann problem consists of a mirror image with respect to thesolid-wall. In the solid-wall boundary conditions, where subscript Lmeans the quantities ρ, p, u, v, θ in the ghost cells; while subscript Rmeans any ones in the boundary cells. θ_(w) is the angle of thesolid-wall. Based on the solutions of the Riemann problem given in(37)-(40), on the solid-wall,P _(*) =P _(R),  (58)

$\begin{matrix}{\rho_{*L} = {\rho_{*R} = {\rho_{R} + {\frac{1}{2}\left( {{f\left( {V_{R},\theta_{L}} \right)} - {f\left( {V_{R},\theta_{R}} \right)}} \right){\left( \frac{\rho_{R}}{a_{R}} \right).}}}}} & (59)\end{matrix}$

The pressure in the star state can be considered as the pressure on thesolid-wall p_(w)P _(w) =P _(*).  (60)

If the pressure on the stationary wall (p_(w)) had been specified, thespecified pressure can be directly used as the pressure in the starstate (region) in the Riemann solver,P _(*) =P _(w).  (61)

In/Out Flow Boundary Condition

For the subsonic inlet boundary, the first relation corresponds to thepositive, incoming, characteristics. Hence the values at the boundary,indicated by a subscribe b, are obtained from the Riemann invariance in(21)-(22). They are

b - = u b + V b V b 2 + U B 2 ⁢ 2 ⁢ a b γ - 1 = u ∞ + V ∞ V ∞ 2 + U ∞ 2 ⁢ 2⁢a ∞ γ - 1 , ( 62 ) b + = u b - V b V b 2 + U B 2 ⁢ 2 ⁢ a b γ - 1 = u i - Vi V i 2 + U i 2 ⁢ 2 ⁢ a i γ - 1 , ( 63 )where u_(∞), a_(∞) are the velocity and the speed of sound for the freestream; the subscript i refers to a value at an internal point. Thesecond relation is associated with a numerical boundary condition andhas to be estimated from inside the domain by an appropriateextrapolation. Hence, the velocity u_(b) is obtained by adding equation(62) and (63), then

u b = b + + i - 2 . ( 64 )

A quite similar way can be used to treat the subsonic outflowboundaries, with the difference that the one specified variable at theboundary is p_(∞), (the pressure at infinite),

b - = -  u b  + V b V b 2 + U B 2 ⁢ 2 ⁢ a b γ - 1 = ∞ - = -  u ∞  + V∞ V ∞ 2 + U ∞ 2 ⁢ p ∞ ρ ∞ ⁢ a ∞ . ( 65 )

Boundary Conditions for Lagrangian Geometry State Variables

Another boundary conditions need to implement are the Lagrangiangeometry state variables U and V. For the solid-wall boundary,U _(L)=sin(2θ_(w)−θ_(R)),  (66)V _(L)=cos(2θ_(w)−θ_(R)).  (67)

For in/outlet boundaries, we specify the flow angle θ=0, then,U=0,  (68)V=1.  (69)

A flow chart presenting the numerical procedures for one step ofupdating from [Q]_(i,j) ^(n) to [Q]_(i,j) ^(n+1) with second-orderaccuracy in time and space, where Q=[ρ, u, v, p]^(T) are the physicalparameters of the flow field, is given in FIG. 2.

Method for Solving a Class of Inverse Problems

Since using the Lagrangian formulation, the solid-wall boundaries arealways represented by streamlines, which will be the stream function inLagrangian plane, where any streamlines represented by the computinggrids are straight lines. Exploring the advantage of the Lagrangianformulation on a class of inverse problem is another purpose in thisinvention. From now, this class problem is defined as the inversegeometry shape design problem.

In Section 2, we introduced the solutions of the local Riemann problemin ξ-direction in Lagrangian plane and the solid-wall boundarycondition. From (37) and (57), ones have

$\begin{matrix}{{P_{w} = {P_{R} + {\frac{\rho_{R}a_{R}}{2}\left( {{f\left( {V_{R},\theta_{L}} \right)} - {f\left( {V_{R},\theta_{R}} \right)}} \right)}}},} & (70)\end{matrix}$where p_(R), ρ_(R), a_(R), V_(R), θ_(R) are the known state variables inboundary and f is the combination function in triangle form given in(35) and (36). From (70), ones have

$\begin{matrix}{{{f\left( {V_{R},\theta_{L}} \right)} = {{f\left( {V_{R},\theta_{R}} \right)} + \frac{2\left( {P_{w} - P_{R}} \right)}{\rho_{R}a_{R}}}},} & (71)\end{matrix}$which can be numerically solved with the specified p_(w) to get the flowangle θ_(L) in ghost cells. Using the mirror boundary condition, thewall angle is

$\begin{matrix}{\theta_{w} = {\frac{1}{2}{\left( {\theta_{L} + \theta_{R}} \right).}}} & (72)\end{matrix}$

The inverse geometry shape design means finding the solid-wall anglebased on the specified pressure distribution on the solid-wall. Therecovering of the solid-wall angle needs to solve the inversed Riemannproblem here. An inverse Riemann solver is built as the followingprocedure. One can definetgθ _(L)=ε,  (73)

Then we have the velocity components in ghost cells

$\begin{matrix}{{u_{L} = {{V_{R}\cos\;\theta_{L}} = \frac{V_{R}}{\sqrt{1 + ɛ^{2}}}}},{v_{L} = {{V_{R}\sin\;\theta_{L}} = \frac{ɛ\; V_{R}}{\sqrt{1 + ɛ^{2}}}}},} & (74)\end{matrix}$which can be substituted into the combination function in triangle form(71), which can be further modified and finalized as a function of ε

$\begin{matrix}{{F(ɛ)} = {{\frac{V_{R}}{2}\left( {ɛ + \frac{{\ln\; V_{R}} + {\ln\left( {ɛ + \sqrt{1 + ɛ^{2}}} \right)} - {\ln\sqrt{1 + ɛ^{2}}}}{\sqrt{1 + ɛ^{2}}}} \right)} - {f\left( {V_{R},\theta_{R}} \right)} - {\frac{2\left( {P_{w} - P_{R}} \right)}{\rho_{R}a_{R}}.}}} & (75)\end{matrix}$

The work needs to do is numerically to solve the equation F(ε)=0 with agiven physical parameters ρ_(R), p_(R), V_(R), θ_(R) and p_(w) to find εfor θ_(L). The equation (75) does have the differentiability and itsfirst-order derivative of F in respect to ε is

$\begin{matrix}{{F^{\prime}(ɛ)} = {\quad{\frac{V_{R}}{2}{\left( \frac{\begin{matrix}{{\left( {2 + ɛ^{2}} \right)\sqrt{1 + ɛ^{2}}} - ɛ -} \\{ɛ\left( {{\ln\; V_{R}} + {\ln\left( {ɛ + \sqrt{1 + ɛ^{2}}} \right)} - {\ln\sqrt{1 + ɛ^{2}}}} \right)}\end{matrix}}{\sqrt{\left( {1 + ɛ^{2}} \right)^{3}}} \right).}}}} & (76)\end{matrix}$

Similar to the method for (45), the updated value ε in new step can beexpressed as

$\begin{matrix}{ɛ^{n + 1} = {ɛ^{n} - {\frac{F\left( ɛ^{n} \right)}{F^{\prime}\left( ɛ^{n} \right)}.}}} & (77)\end{matrix}$

The flow angle in ghost cells can be recovered by θ_(L)=tg⁻¹(ε^(n+1)).The solid-wall angle θ_(w) can be found by using equation (72).

The flow chart for solving the inverse geometry design problem is givenin FIG. 3.

This method starts from the guessed flow filed physical values andsolid-wall angle, the inversed Riemann solver with the specifiedpressure distribution applied on the solid wall boundary will give theflow angles in ghost cells, then the wall angle can be obtained andsubstituted back to the flow solver for updating the physical variables.This procedure will continue until an iterated solid-wall angle withinthe convergence tolerance is reached. It is obvious that the convergentcriteria is the solid-wall angle instead of any physical parameter,which means that the solid-wall angle become a flow field variable. Thechange of the boundary each iteration doesn't change the characteristicsof the equations and the equations are still the strong-hyperbolicsystem; therefore, it is still a well-posed initial problem. Compared tothe methods used in the Eulerian plane, such as the adjoint method, themethod needs neither the step of the grid generation for the updatedgeometry shape, nor the step of solving the adjoint equation. It obtainsthe convergent physical solutions in flow field and the geometry shapeconcurrently, which will be the one-shot and optimal method in theinverse geometry shape design problem. The total CPU time for this classproblem will reduce the order of one.

Invicid Subsonic Internal Flow Passing a Divergent Nozzle

According to the numerical methods discussed above, a complete solutionto the invicid subsonic internal flow by using the Euler equations inLagrangian formulation will be performed. This example is about thesubsonic flow passing through a two-dimensional parabolic divergentnozzle with length L and inlet height H_(in), whose geometry is definedas the two parts of parabolic curves

$\begin{matrix}{{H(x)} = \left\{ \begin{matrix}{{- {ax}^{2}},} & {{0 \leq x < \frac{L}{2}};} \\{{{a\left( {x - L} \right)}^{2} - b},} & {{\frac{L}{2} \leq x \leq L},}\end{matrix} \right.} & (78)\end{matrix}$where H_(in)=L/3, a=H_(in)/2, b=H_(in)/L. The Mach number of theinviscid flow at the nozzle inlet is M_(in)=0.5. The flow is the puresubsonic flow. The Euler equations in Lagrangian formulation (7) and (8)will be solved by the numerical scheme presented in Section 2, where theStrang dimensional-splitting of second-order accuracy in time and hybridflux operators with second-order upwind MUSCL extrapolation with minmodlimiter are used. The computing parameters: CFL=0.48; the Lagrangianbehavior parameter h=0.98; totally 60×20 uniform cells in Lagrangianplane. The computing results will be compared with those obtained fromthe JST method based on FVM in Eulerian plane and the exact solutions.Now, denoting the solution obtained in Lagrangian plane by the inventedmethod as the Lagrangian solutions and the results obtained from theEulerian plane by the Eulerian formulation as the Eulerian solutions. Indetail, the computations are taken as the following steps with referringto FIG. 2.

(1) Initialization. The flow field variables Q⁰=[ρ⁰, u⁰, v⁰, p⁰]^(T) areknown as the initial conditions as well as U⁰, V⁰, where Q⁰ take thevalues at the infinite and U⁰=0; V⁰=1. The superscript 0 means that theinitial conditions of a flow problem are given at t=0 (τ=0) in x-y planeand λ-ξ plane. Subsequently, the conservation variables f⁰ areavailable. Then an appropriate rectangular λ-ξ coordinate mesh withuniform spatial-step (Δλ, Δξ) is formed. The corresponding mesh in x-yplane is laid on according todx=cos θdλ and dy=sin θdλ.  (79)

(2) Calculating the time-step. Based on the variables at time-step n,(n=0, 1, 2, . . . ), and CFL<0.5, calculate the time-step,

$\begin{matrix}{{{\Delta\;\tau} = {{CFL} \cdot {\max\left( {\frac{\Delta\;\lambda}{{{\left( \frac{a}{J} \right)\sqrt{U^{2} + V^{2}}} + {\left( {1 - h} \right)\sqrt{u^{2} + v^{2}}}}},\frac{\Delta\xi}{\frac{a}{J}}} \right)}}},} & (80)\end{matrix}$where S=sin θ^(n−1), T=cos θ^(n−1), where θ^(n−1) is from previoustime-step, thenj=UT−VS.  (81)

(3) Using the FVS method in λ-direction to solve equations (23). Theconservation variables are updated from f_(i,j) ^(n) to f_(i,j) ^(n+1/2)by the known values [ρ, u, v, p]_(i,j) ^(n) and U_(i,j) ^(n), V_(i,j)^(n). For the second-order schemes, the MUSCL extrapolation with TVDlimiter has to be taken.

(4) Using the Godunov method in ξ-direction to solve equation (29). Withall f^(n+1/2) and Q^(n+1/2) known at time step n (n=0, 1, 2, . . . ),the Godunov method will be used to update the conservation variables. Acomplete 2-D local Riemann solver gives the primitive variables [ρ, u,v, p]_(i,j±1/2) at the interfaces of the neighboring cells. For thesecond-order schemes, the MUSCL extrapolation with TVD limiter has to betaken to decide the initial values for the Riemann problem.

(5) Updating the geometry variables. Since the velocities [u,v]_(i,j±1/2) ^(n+1/2) at the interface of cells are available, then thediscretized geometrical conservation equations can be updated to newtime step at once by (57), whose accuracy in time and space is decidedby the accuracy of the velocities [u, v]_(i,j±1/2) ^(n+1/2) though thediscretized equations are presented in the first-order accuracy.

(6) Finding the physical primitive variable values at the new time step.Decoding from the conservation variables f_(i,j)=[f₁, f₂, f₃, f₄]_(i,j)^(n+1) to obtain the primitive variables [ρ, p, u, v]_(i,j) ^(n+1). Theprimitive variables read as

$\begin{matrix}{{\rho = \frac{f_{1}}{J}},{u = \frac{f_{2}}{f_{1}}},{v = \frac{f_{3}}{f_{1}}},{p = {\frac{\gamma - 1}{J}{\left( {f_{4} - \frac{f_{2}^{2} + f_{3}^{2}}{2f_{1}}} \right).}}}} & (82)\end{matrix}$

(7) Grid building. The last step in one time-step is to build the gridif needed. The grid points are given byx _(i,j) ^(n+1) =x _(i,j) ^(n)+½Δτ_(i,j) ^(n)(hu _(i,j) ^(n) +hu _(i,j)^(n+1)); y _(i,j) ^(n+1) =y _(i,j) ^(n)+½Δτ_(i,j) ^(n)(hv _(i,j) ^(n)+hv _(i,j) ^(n+1))  (83)according to which each grid point represents the center of a fluidparticle representing by a computational cell. The principle of thetime-dependent Lagrangian approach in this invention tells that alongξ-direction, the grid lines built here actually are the streamlinethemselves.

(8) Loop (2). After the numerical procedure for advancing one time-stepis complete, to march forward further, one goes back to (2) and repeat(3)-(7), until the conservation variables or primitive variables go toconvergence.

The Lagrangian solution starts from the rectangular grid with Lagrangiandistance in λ-direction and stream function in ξ-direction in FIG. 4(a). The streamlines shown in FIG. 4( c) by Lagrangian solutionapparently agreed well with Eulerian solution in FIG. 4( d), and thepressure distributions on the solid-wall shown in FIG. 4( e) totallyagreed with the exact solution. The grid used in the Euleriancoordinates in FIG. 4( b) has little difference with the grid built bythe Lagrangian solution, where the lines along x-direction arestreamlines and the lines along y-direction are not perpendicular tox-direction to preserve the length of the streamlines. It is worth toindicate that the final grids in FIG. 4( f) evolve from the initial gridin FIG. 4( a), and it is the result of the Lagrangian solution since thegrid lines are the streamlines themselves. Although one does not needthis grid in the computing, an evolving procedure can give a clearunderstanding of the evolution of the Lagrangian solution.

Inverse Geometry Shape Design Case

Next case is an example to show the detail how to find the geometry of asolid-wall by solving the two-dimensional Euler equation in Lagrangianformulation. It is an inverse shape design problem. The geometry of theexample is a convergent-divergent nozzle. Since the prescribed pressuredistribution on the geometrically-unspecified solid-wall need to beinput in the inverse problem, one can firstly specify this nozzle'sshape and find the pressure distribution, then use this pressure as theinput to find the geometry shape by the invented method. The calculatedshape should match with the original shape. The profile of the nozzle'ssolid-wall is given by a functiony=H _(th) e ^(−βx) ² ,−1≦x≦1,  (84)where H_(th) is taken to 0.1 meaning the 10% bump high of the inlet andβ is taken 6.5 to control the smoothed convex part which counts about athird of the channel length. For an inlet Mach number of 0.5, it is apure subsonic flow. The grid generation of this geometry in theLagrangian plane is identical to that in the previous case shown in FIG.4( a). The computational grid in Euler plane is given in FIG. 5( a). TheEuler equations in Lagrangian formulation (7) and (8) will be solved bythe numerical scheme presented in Section 2, where the Strangdimensional-splitting of second-order accuracy in time and hybrid fluxoperators with second-order upwind MUSCL extrapolation with minmodlimiter are used. The computing parameters: CFL=0.48; the Lagrangianbehavior parameter h=0.98; and two kinds grid, totally 60×20 and 90×30uniform cells in Lagrangian plane. To obtain the specified pressuredistributions on the solid-wall, the flow field calculations in theLagrangian and Eulerian plane firstly are implemented like in theprevious case. FIG. 5( b) shows the pressure distribution on thesolid-wall by the Lagrangian and the Eulerian solution, which will beused as the prescribed pressure distribution in the inverse geometryshape design problem. Also, the adjoint method will be used to solvethis inverse design problem to test the computing efficiency of theinvented method. The computations are taken as the following steps, indetail:

(1) Initializing the flow field and the solid-wall angle. The initialedflow field variables are set as the same as in the previous case. Thesolid-wall angle is assumed as θ_(w)=0⁰. Since the flow angle on thewhole flow field are assumed as zero, the flow angle in the ghost cellsis zero θ_(w)=0⁰.

(2) Imputing the prescribed pressure distribution on the solid-wallp_(w).

(3) Solving the inversed Riemann problem across the streamline on thesolid-wall boundary based on p_(w), according to the equations (70)-(77)in Section 2.

(4) Calculating the solid-wall angle. The solid-wall angle can becalculate based on the mirror state condition, then we haveθ_(w)=½θ_(L)+θ_(R)).  (85)

$\begin{matrix}{\theta_{w} = {\frac{1}{2}{\left( {\theta_{L} + \theta_{R}} \right).}}} & (85)\end{matrix}$

(5) Checking the convergence of the solid-wall angle. The variations ofthe solid-wall angle is the difference of θ_(w) between the time stepn+1 and nδθ_(w)=(θ_(w) ^(n+1)−θ_(w) ^(n))  (86)

If δθ_(w) is smaller than the convergent tolerance, then the flow fieldand solid-wall angle come to the final solution θ_(w) ^(n+1). Otherwise,continue the next steps.

(6) Updating the flow field. The task of this step is to calculate theflow field physical parameters from [Q]_(i,j) ^(n) to [Q]_(i,j) ^(n+1)based the updated solid-wall angle θ_(w) ^(n+1). This step includes thestep (1)-(7) introduced in the previous case. After this step, the newflow field parameters (ρ_(R), p_(R), V_(R), θ_(R)) are obtained.

(7) Repeating (2).

FIG. 5( c) shows the graphic comparison of the computed and originalsolid-wall shapes. From the results, we can numerically evaluate theprecision of the resolution. The computed geometrical shape was found ingood agreement with the original shape used as the input; the maximumerror is limited within 0.5%. The totally CPU time is only tenth of thatby the adjoint method.

While particular embodiments of the invention have been shown anddescribed, it will be obvious to those skilled in the art that changesand modifications may be made without departing from the invention inits broader aspects, and therefore, the aim in the appended claims is tocover all such changes and modifications that fall within the truespirit and scope of the invention.

REFERENCE

-   (1) Godunov, S. K.: A difference scheme for numerical computation    discontinuous solution of hydrodynamic equations. Translated US    Publ. Res. service, JPRS 7226, 1969.-   (2) Van Leer, B.: Flux vector splitting for the 1990's. Invited    Lecture for the CFD Symposium on Aero propulsion, Cleveland, Ohio,    1990.-   (3) Jameson, A.: Aerodynamic design via control theory. Journal of    Scientific Computing, 3:233-260, 1988.15.-   (4) Loh, C. Y. and Liu, M. S.: A new Lagrangian method for    three-dimensional steady supersonic flows. Journal of Computational    Physics, 113, pp. 224-248, 1994.-   (5) Mateescu, D.: Analysis of aerodynamic problem with geometrically    unspecified boundaries using an enhanced Lagrangian method. Journal    of Fluids and Structures, 17, pp. 603-626, 2003.-   (6) Stanitz, D. John: A review of certain inverse methods for the    design of duct with 2- or 3-dimensional potential flow. Appl. Mech.    Rev. Vol. 41, No. 6, June 1988.-   (7) Jameson, A., Schmidt, W., Turkel, E.: Numerical solutions of the    Euler equations by finite volume methods using Runge-Kutta    time-stepping schemes. AIAA-paper, 81-1259, 1981.

The invention claimed is:
 1. A computer implemented method forsimulating two-dimensional inviscid subsonic flows, comprising thefollowing steps: (1) Transforming, by using a computer, a set oftwo-dimensional unsteady Euler equations in a Euler plane into a set ofunsteady Euler equations in Lagrangian formulation in a Lagrangianplane, which has one time coordinate direction stated by a Lagrangiantime τ and two space coordinate directions stated by a stream function ξand a Lagrangian-distance λ, through a transforming matrix${J = \begin{bmatrix}1 & 0 & 0 \\{hu} & {\cos\;\theta} & U \\{hv} & {\sin\;\theta} & V\end{bmatrix}};$  said unsteady Euler equations in Lagrangianformulation include two set equations: a set of Lagrangian physicalconservation equations and a set of geometric conservation equations,which formally are written as: the Lagrangian physical conservationequations${{\frac{\partial f_{L}}{\partial\tau} + \frac{\partial F_{L}}{\partial\lambda} + \frac{\partial G_{L}}{\partial\xi}} = 0},$ where f_(L) is the conservation variables vector and F_(L), G_(L) arethe convective fluxes respectively in the λ and ξ direction in theLagrangian plane, in detailed ${f_{L} = \begin{bmatrix}{\rho\; J} \\{\rho\;{Ju}} \\{\rho\;{Jv}} \\{\rho\;{JE}}\end{bmatrix}},{F_{L} = \begin{bmatrix}{\left( {1 - h} \right)\rho\; K} \\{{\left( {1 - h} \right)\rho\;{Ku}} + {Vp}} \\{{\left( {1 - h} \right)\rho\;{Kv}} - {Up}} \\{\left( {1 - h} \right)\rho\;{KH}}\end{bmatrix}},$ ${G_{L} = \begin{bmatrix}0 \\{{- p}\;\sin\;\theta} \\{p\;\cos\;\theta} \\0\end{bmatrix}},{\theta = {{tg}^{- 1}\left( \frac{v}{u} \right)}},$K=uV −vU, 0<h≦1; the geometric conservation equations${{\frac{\partial f_{g}}{\partial\tau} + \frac{\partial F_{g}}{\partial\lambda} + \frac{\partial G_{g}}{\partial\xi}} = 0},{{{where}\mspace{14mu} f_{g}} = \begin{bmatrix}U \\V\end{bmatrix}},{F_{g} = \begin{bmatrix}0 \\0\end{bmatrix}},{G_{g} = {\begin{bmatrix}{- {hu}} \\{- {hv}}\end{bmatrix}.}}$  in all above, the variables t, ρ, p, E and H arerespectively the time, fluid density, pressure, total energy andenthalpy; u, v are the Cartesian components of flow velocity V, and U, Vare the Lagrangian geometry state variables; (2) reading object geometryfor providing points on a surface of an object; (3) establishing acomputing mesh around said object; (4) using aStrang-dimensional-splitting scheme with hybrid upwind flux operatorssolving the unsteady Euler equations in Lagrangian formulation in theLagrangian plane.
 2. The method of claim 1, wherein said computing meshis a two-dimensional rectangular mesh with λ and ξ as two coordinatedirections in the Lagrangian plane.
 3. The method of claim 1, whereinsaid using the Strang-dimensional-splitting scheme with hybrid upwindflux operators solving the unsteady Euler equations in Lagrangianformulation in the Lagrangian plane comprises the following steps in onetime step for updating the conservative variables f_(L) when theLagrangian time τ in the Lagrangian plane evolves, (1) calculating theconvective fluxes across cell interfaces of the computing mesh inλ-direction by the FVS method with a time-step Δτ; (2) solving a Riemannproblem across streamlines in ξ-direction with half time-step Δτ by theGodunov method with a local Riemann solver sequentially; (3) calculatingthe convective fluxes across the cell interfaces in λ-direction by theFVS method with the time-step Δτ again.
 4. The method of claim 3,wherein said Riemann problem across streamlines in ξ-direction comprisesthat at a left or right state (region) of the cell interface, are shockand expansion waves respectively, between which it is a star state(region), which is divided into a left and a right star states (region)again.
 5. The method of claim 3, wherein said solving the Riemannproblem across streamlines comprises the following steps: (1) Connectinga left and a right states to the star state by integrating along a setof characteristic equations of the Lagrangian physical conservationequations; (2) Recovering velocity magnitude in the star state; (3)Solving a combination function to find a flow angle in the star state;(4) Finding the velocity components in the star state.
 6. The method ofclaim 5, wherein said step of recovering the velocity magnitude in thestar region is done according to a Rankine-Hugoniot relation across theshocks and an Enthalpy constants relation across the expansion waves. 7.The method of claim 5, wherein said combination function f(u, v) iswritten as${f\left( {u,v} \right)} = {{\frac{1}{2}\frac{v\sqrt{u^{2} + v^{2}}}{u}} + {\frac{1}{2}u\;{{\ln\left( {v + \sqrt{u^{2} + v^{2}}} \right)}.}}}$